Trabajos Prácticos de Simulación [75.26]¶

  • Albarracin, Dario Javier - 97568
  • Bergman, Guido Ernesto - 104030
  • del Mazo, Federico - 100029
  • Rodriguez, Florencia - 100033
  • Salgueiro Da Costa Riberio, Joao Luís

Trabajo Práctico 1¶

LXM¶

LXM: Better Splittable Pseudorandom Number Generators (and Almost as Fast)

  • Paper: https://dl.acm.org/doi/abs/10.1145/3485525
  • Video: https://youtu.be/OXurCqln_qc

En el 2021 Guy Steele y Sebastiano Vigna presentaron un generador de números pseudo aleatorios (PRNG) que se construye a partir de otro PRNG en donde colaboró Guy Steele llamado SplitMix.

La principal propiedad de ambos PRNGs es que estos generadores son partibles (splittable), es decir, se pueden dividir en dos nuevos generadores estadísticamente independientes, lo cual es de gran utilidad en ambientes concurrentes.

LXM se basa en combinar tres ideas: un generador lineal congruente (la L), un generador basado en XORs (la X) y el resultado de la combinación de esos dos generadores utilizarlo como input en una función mezcladora (la M).

En este trabajo nos centraremos específicamente en la parte práctica del algoritmo y lo re-implementaremos en nuestro propio código, evitando hablar tanto de la concurrencia como de la propiedad partible del algoritmo (lo cual consiste de la mayoría del paper de Guy Steele y Sebastiano Vigna).

Las ideas principales de este trabajo surgen de la sección 2 del paper (THE LXM GENERATION ALGORITHM)

LXM.png

In [1]:
def LXM(s, m=2891336453, a=1310709051, x0=11565345, x1=5242836, k=64):
    
    # https://gist.github.com/trietptm/5cd60ed6add5adad6a34098ce255949a
    rotl = lambda val, r_bits, max_bits: \
        (val << r_bits%max_bits) & (2**max_bits-1) | \
        ((val & (2**max_bits-1)) >> (max_bits-(r_bits%max_bits)))    
    
    while True:
        # Combining operation
        z = s + x0
        # Mixing function (lea64)
        z = (z ^ (z >> 32)) * 0xdaba0b6eb09322e3
        z = (z ^ (z >> 32)) * 0xdaba0b6eb09322e3
        z = (z ^ (z >> 32))
        z &= (2**k) - 1
        # Update the LCG subgenerator
        s = (m * s + a)
        # Update the XBG subgenerator (xoroshiro128v1_0)
        q0 = x0; q1 = x1
        q1 ^= q0
        q0 = rotl(q0, 24, k)
        q0 = q0 ^ q1 ^ (q1 << 16)
        q1 = rotl(q1, 37, k)
        x0 = q0; x1 = q1
        # Return result (modified to be in interval 0-1)
        yield z/2**k
In [2]:
n = 5
print(f"{n} números generados al azar, en el intervalo 0-1:")

g = LXM(42) 
for r, i in zip(g,range(n)):
    print(f"* {r}")
5 números generados al azar, en el intervalo 0-1:
* 0.793781427453092
* 0.2155520740899848
* 0.7433470653358116
* 0.28660629734277526
* 0.515362051632894
In [5]:
scatterplot.show()

En el gráfico de dispersión se pretende ver si la generación de los números sigue algún patrón. En el gráfico realizado no se observa ninguno, ya que no se ven formas claras, como por ejemplo una línea recta, ni zonas con una densidad marcadamente superior a la del resto. Lo que se observa es una "nube de puntos", que es lo que debería verse en una distribución uniforme.

In [7]:
scatter3D.show()
Text(0.5, 0, 'X(n)+2')
In [9]:
histogramas.show()

En un histograma se muestra la cantidad de observaciones de la variable estudiada según el intervalo. Permite observar como se distribuyen los números generados. Por ejemplo, si las barras graficadas tuvieran diferencias significativas en sus alturas esto querría decir que la distribución no es uniforme. En los histogramas hechos, podemos ver que a medida que aumenta la cantidad de muestras usada, y consecuentemente se reduce el tamaño del intervalo, la diferencia que se observa entre las barras disminuye. Esto culmina en el último gráfico, que usa 52016 muestras e intervalos de tamaño 0.06, en el que la longitud de las barras es prácticamente igual. Por esta razón podemos decir que utilizando una cantidad de muestras suficientes (ej: 52016), las muestras tenderán a distribuirse uniformemente en los intervalos.

In [11]:
boxplot.show()

Un boxplot muestra los cuartiles de los números generados, el soporte de la distribución y los outliers. Para el caso actual, se puede utilizar para comparar esos datos con los de una distribución uniforme. Así, se puede ver que:

  • Los cuartiles son los esperados. Esto se debe a que el del 25% coincide con 0.25, el de 50% con 0.5 y el de 75% con 0.75
  • El soporte es entre 0 y 1
  • No se observan outliers
In [13]:
qqplot.show()

En el QQplot se comparan los cuantiles de la variable a estudiar, en este caso los de los números generados con LXM, con los de una variable ideal, en este caso una distribución uniforme. Este gráfico permite observar si la probabilidad acumulada a izquierda de cada valor de la variable a estudiar se aleja o no del de la variable ideal. Se puede ver que si bien en el primer gráfico se observan claras diferencias entre los cuantiles teóricos y los de la muestra, estas pueden ser atribuidas a la baja cantidad de muestras tomadas para dicho gráfico. Lo anterior se debe a que dichas diferencias desaparecen al aumentar la cantidad de muestras en los gráficos siguientes. De hecho, se puede apreciar que en el gráfico hecho utilizando 1000 muestras, aproximadamente entre 0.1 y 0.4 los cuantiles de la muestra se encuentran apenas por debajo de los teóricos, pero esto se corrige en el gráfico hecho 10000 muestras. Por lo tanto, podemos decir que tomando una cantidad de muestras lo suficientemente grande (ej:10000), los cuantiles de los números generados y los de una distribución aleatoria no presentan diferencias.

Test Estadísticos¶

  • Test de Chi2
  • Test de Gap
  • Test de Independencia

Test de Chi2¶

Con este método comparamos dos distribuciones de probabilidad. Por una parte, usaremos el generador que implementamos, y por otra parte, para la distribución esperada.

Con este método comparamos dos distribuciones de probabilidad. Por una parte, la observada, que es la que obtenemos utilizando nuestro generador LXM. Por otra parte, la distribución esperada, que es la distribución que queremos corroborar que la observada sigue.

In [14]:
from scipy.stats import chi2
import numpy as np
import math

class Chi2:
    def __init__(self, numero_de_muestras=1000, g = None):
        self._generador = g or LXM(42)
        self._numero_de_muestras = numero_de_muestras
        
    def realizar_prueba(self):
        frecuencias_observadas = self.frecuencias_observadas()
        
        f_esperada = 1/10 * sum(frecuencias_observadas)
        d2 = sum([(f_observada - f_esperada)**2 for f_observada in frecuencias_observadas]) / f_esperada

        limite_superior = chi2.ppf(0.95, df=5)
        print(f'Frecuencias: {frecuencias_observadas}')
        print('Estadístico: {:.2f}'.format(d2))
        print(f'Límite superior: {round(limite_superior, 2)}')
        self.validar_hipotesis(d2, limite_superior)

    def calcular_d2(self):
        numerador = (f_observada - f_esperada) ** 2
        denominador =  f_esperada
        d_al_cuadrado = numerador / denominados # acá es sumatoria
    
    def validar_hipotesis(self, d2, limite_superior):
        if d2 <= limite_superior:
            print("El test acepta la hipotesis nula")
        else:
            print("El test rechaza la hipótesis nula")
        
    def frecuencias_observadas(self):
        # Estas son las que tengo que obtener
        numeros_generados = []
        for i in range(self._numero_de_muestras):
            numeros_generados.append(next(self._generador))
        frecuencias = [0] * 10
        for numero in numeros_generados:
            pos = math.floor(numero*10)
            if pos < 10:
                frecuencias[pos] += 1
        return frecuencias
In [15]:
def pseudo_chi2():
    # Obtengo freuencias observadas de 0 a 9 (enteros)
    # La frecuencia esperada es 1/10  por la cantidad de números aleatorios que genero
    # Con estos valores ahora entonces calculo d2
    # d2 = sum (frecuencia_observada - frecuencia_esperada) ** 2  / frecuencia_esperada
    # Finalmente para chi2.ppf(0.95, df=5) vemos si d2 es mayor o menos al limite superior 
    return
In [16]:
# Corremos nuestro experimento
Chi2(numero_de_muestras=10000).realizar_prueba()
Frecuencias: [949, 998, 1019, 1007, 1008, 998, 1029, 968, 1001, 1023]
Estadístico: 5.48
Límite superior: 11.07
El test acepta la hipotesis nula

Test de Gap¶

Este test consiste en contar la cantidad de números aleatorios generados que no pertenecen a un intervalo dado, tal que consigamos 'gaps'. Para este test decidimos generar 1000 gaps y utilizar el intervalo [0.2, 0,5].

Finalmente aplicamos chi2 para validar (o rechazar) la hipótesis.

In [17]:
import matplotlib.pyplot as plt

class Gap:
    def __init__(self, cota_inferior, cota_superior, g = None):
        self._generador = g or LXM(42)
        self._cota_inferior = cota_inferior
        self._cota_superior = cota_superior
        self._p_0 = cota_superior - cota_inferior
        self._cantidad_de_gaps = 1000
        
    def realizar_prueba(self):
        # Espero a la primera vez que tengo valor dentro de la cota
        gaps = self.generar_gaps()
        self.plotear_gaps(gaps)
        frecuencias_observadas = self.generar_frecuencias(gaps)
        d2 = 0
        
        # Aplicamos chi2
        for i, frecuencia_observada in enumerate(frecuencias_observadas):
            f_esperada = self.frecuencia_esperada(i)
            numerador = (frecuencia_observada - f_esperada) ** 2
            denominador = f_esperada
            d2 += numerador / denominador
        
        # Validamos hipótesis
        k = len(frecuencias_observadas)
        limite_superior = chi2.ppf(0.95, df=k-1)
        print('Estadístico: {:.2f}'.format(d2))
        print(f'Límite superior: {round(limite_superior, 2)}')
        self.validar_hipotesis(d2, limite_superior)
        
    def validar_hipotesis(self, d2, limite_superior):
        if d2 <= limite_superior:
            print("El test acepta la hipotesis nula")
        else:
            print("El test rechaza la hipótesis nula")
            
    def generar_gaps(self):
        gaps = []
        # Generamos 1000 gaps
        while(len(gaps) < self._cantidad_de_gaps):
            while (True):
                n = next(self._generador)
                if n < self._cota_superior and n > self._cota_inferior:
                    break
            gap = 0
            while(True):
                n = next(self._generador)
                if n < self._cota_superior and n > self._cota_inferior:
                    gaps.append(gap)
                    break
                else:
                    gap +=1
        return gaps
    
    def plotear_gaps(self, gaps):
        plt.hist(gaps, color='b', bins="sturges", alpha=0.8, density=True, rwidth=0.8, align='left')
        plt.show()
        
    def generar_frecuencias(self, gaps):
        frecuencias = [0] * (max(gaps) + 1)
        for gap in gaps:
            frecuencias[gap] += 1
        return frecuencias
    
    def frecuencia_esperada(self, n):
        return (((1 - self._p_0)**n) * self._p_0) * self._cantidad_de_gaps
        
In [18]:
def pseudo_gap():
    # Obtengo gaps entre un rango de valores 
    # (dentro de [0,1]  ya que son los valores que genera el LXM programado
    # Con estos valores ahora entonces calculo d2
    # d2 = sum (frecuencia_observada - frecuencia_esperada) ** 2  / frecuencia_esperada
    # La  frecuencia esperada podrá ser distinta para cada valor, ya que
    # la misma sigue una distribución geométrica
    # Finalmente para chi2.ppf(0.95, df=cantidad_de_gaps) vemos si d2 es mayor o menos al limite superior 
    return
In [19]:
# Corremos nuestro experimento
Gap(cota_inferior=0.2, cota_superior=0.5).realizar_prueba()
Estadístico: 21.92
Límite superior: 31.41
El test acepta la hipotesis nula

Test de Independencia¶

Con este métod vamos a validar que las variables sean independientes.

In [20]:
from scipy.stats import chi2_contingency


class TestIndependencia:
    def __init__(self, g = None):
        self._generador = g or LXM(42)
        self._cantidad_de_vectores = 500
        
    def realizar_prueba(self):
        nros_x = []
        nros_y = []
        
        for i in range(self._cantidad_de_vectores):
            nros_x.append(next(self._generador))
            nros_y.append(next(self._generador))
        
        self.plot_vectores(nros_x, nros_y)
        
        d2, p, dof, ex = chi2_contingency([nros_x, nros_y], correction=True)
        
        # Validamos hipótesis
        limite_superior = chi2.ppf(0.95, df=dof)
        print('Estadístico: {:.2f}'.format(d2))
        print(f'Límite superior: {round(limite_superior, 2)}')
        self.validar_hipotesis(d2, limite_superior)
        
    def validar_hipotesis(self, d2, limite_superior):
        if d2 <= limite_superior:
            print("El test acepta la hipotesis nula")
        else:
            print("El test rechaza la hipótesis nula")

    def plot_vectores(self, nros_x, nros_y):
        plt.title(f"Dispersión de {self._cantidad_de_vectores} números generados al azar")
        plt.ylabel("X(n)-1")
        plt.xlabel("X(n)")
        plt.scatter(nros_x,nros_y,50)
        plt.show()
            
In [21]:
# Corremos nuestro experimento
TestIndependencia().realizar_prueba()
Estadístico: 87.60
Límite superior: 552.07
El test acepta la hipotesis nula

Ejercicio 3¶

Utilizando el generador implementado en el ejercicio 1:

  • Implementar un método para generar variables aleatorias con distribución normal con media 10 y desvío 2.
  • Graficar la distribución que siguen los números pseudoaleatorios generados.
  • Realizar, al menos 2 tests, de los explicados en la materia, para verificar si los números generados siguen la distribución pedida (evalué los resultados para distintos tamaños de muestra).
In [23]:
gen = LXM(42)

def generar_exponencial(lambda_):           
    u_1 = next(gen)
    return -np.log(u_1)/lambda_

def cociente(t):                            
    return e**((-(t-1)**2)/2)

def transformar_gaussiana(x, mean, std):    
    return x*std + mean

def generar_normal_0_1():
    no_aceptada = True
    while(no_aceptada):
        exp_1 = generar_exponencial(1)
        prob = cociente(exp_1)                
        if(next(gen) < prob):
            no_aceptada = False
            z = exp_1
            if(next(gen) < 0.5):
                z = -z
    return z
In [24]:
m1, m2, m3 = [], [], []     #muestras
for i in range(100):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m1.append(z2)

for i in range(1000):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m2.append(z2)

for i in range(10000):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m3.append(z2)
In [25]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 5))
fig.tight_layout(pad=3)
ax1.hist(m1, 50)
ax1.set_xlabel('Valores generados', fontsize=15)
ax1.set_title("Histograma (n = 100)", fontsize=20)
ax2.hist(m2, 50)
ax2.set_xlabel('Valores generados', fontsize=15)
ax2.set_title("Histograma (n = 1000)", fontsize=20)
ax3.hist(m3, 50)
ax3.set_xlabel('Valores generados', fontsize=15)
ax3.set_title("Histograma (n = 10000)", fontsize=20)

plt.show();

Test de Kolmogorov-Smirnov¶

El test de Kolmogorov-Smirnov se utiliza para decidir si una muestra proviene de una población con una distribución específica.

Se tiene una distribución de probabilidad acumulada, la cual se propone como distribución teórica de la muestra, y otra distribución de probabilidad empirica, generada a partir de la muestra.

En este caso tenemos las hipótesis:



$H_0$: La muestra pertenece a una poblacion con la distribución teórica propuesta
$H_1$: La muestra no pertenece a una poblacion con la distribución teórica propuesta.

El test provee como resultados un estadístico y un p-valor. Cuando el p-valor obtenido es superior a cierto umbral --usualmente 0.05--, decimos que no podemos rechazar la hipótesis de que los datos provienen de una población con la distribución teórica propuesta.
In [26]:
m1_ks, m2_ks, m3_ks = [], [], []     #muestras
for i in range(100):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m1_ks.append(z2)

for i in range(1000):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m2_ks.append(z2)

for i in range(5000):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m3_ks.append(z2)
In [27]:
state = np.random.default_rng()
dist_teorica_1 = sp.stats.norm(loc=10, scale=2).rvs(size=100, random_state=state)
dist_teorica_2 = sp.stats.norm(loc=10, scale=2).rvs(size=1000, random_state=state)
dist_teorica_3 = sp.stats.norm(loc=10, scale=2).rvs(size=5000, random_state=state)
In [28]:
figsts, (ax1sts, ax2sts, ax3sts) = plt.subplots(1, 3, figsize=(10, 3))
figsts.tight_layout(pad=3)

yval1 = np.linspace(0.0, 1.0,num=100)
ax1sts.step(x=np.sort(m1_ks), y=yval1)
ax1sts.step(x=np.sort(dist_teorica_1), y=yval1)
ax1sts.set_xlabel('Valores generados (n=100)', fontsize=8)
ax1sts.set_title("Probabilidad Acumulada", fontsize=10)

yval2 = np.linspace(0.0, 1.0,num=1000)
ax2sts.step(x=np.sort(m2_ks), y=yval2)
ax2sts.step(x=np.sort(dist_teorica_2), y=yval2)
ax2sts.set_xlabel('Valores generados (n=1000)', fontsize=8)
ax2sts.set_title("Probabilidad Acumulada", fontsize=10)

yval3 = np.linspace(0.0, 1.0,num=5000)
ax3sts.step(x=np.sort(m3_ks), y=yval3)
ax3sts.step(x=np.sort(dist_teorica_3), y=yval3)
ax3sts.set_xlabel('Valores generados (n=5000)', fontsize=8)
ax3sts.set_title("Probabilidad Acumulada", fontsize=10)
Out[28]:
Text(0.5, 1.0, 'Probabilidad Acumulada')
In [29]:
figks, (ax1ks, ax2ks, ax3ks) = plt.subplots(1, 3, figsize=(20, 5))
figks.tight_layout(pad=3)
ax1ks.hist(m1_ks, 10)
ax1ks.set_xlabel('Valores generados', fontsize=15)
ax1ks.set_title("Histograma (n = 100)", fontsize=20)
ax2ks.hist(m2_ks, 25)
ax2ks.set_xlabel('Valores generados', fontsize=15)
ax2ks.set_title("Histograma (n = 1000)", fontsize=20)
ax3ks.hist(m3_ks, 50)
ax3ks.set_xlabel('Valores generados', fontsize=15)
ax3ks.set_title("Histograma (n = 5000)", fontsize=20)

plt.show();
In [30]:
from scipy import stats
ks1 = stats.kstest(m1_ks, dist_teorica_1, alternative='two-sided')
ks2 = stats.kstest(m2_ks, dist_teorica_2, alternative='two-sided')
ks3 = stats.kstest(m3_ks, dist_teorica_3, alternative='two-sided')

display ("Valor del estadístico del test Kolmogorov-Smirnov con n = 100: " + str(ks1.statistic))
display ("P-valor del test Kolmogorov-Smirnov con n = 100: " + str(ks1.pvalue))
display ("Valor del estadístico del test Kolmogorov-Smirnov con n = 1000: " + str(ks2.statistic))
display ("P-valor del test Kolmogorov-Smirnov con n = 1000: " + str(ks2.pvalue))
display ("Valor del estadístico del test Kolmogorov-Smirnov con n = 5000: " + str(ks3.statistic))
display ("P-valor del test Kolmogorov-Smirnov con n = 5000: " + str(ks3.pvalue))
'Valor del estadístico del test Kolmogorov-Smirnov con n = 100: 0.07'
'P-valor del test Kolmogorov-Smirnov con n = 100: 0.9684099261397212'
'Valor del estadístico del test Kolmogorov-Smirnov con n = 1000: 0.038'
'P-valor del test Kolmogorov-Smirnov con n = 1000: 0.4659595288557257'
'Valor del estadístico del test Kolmogorov-Smirnov con n = 5000: 0.0178'
'P-valor del test Kolmogorov-Smirnov con n = 5000: 0.40671816589876386'

Dado que en todos los casos el p-valor es superior a 0.05, no rechazamos la hipótesis de que los datos provienen de una distribución normal.

Test de Shapiro-Wilk¶

El test de Shapiro-Wilk está diseñado para decidir si una muestra aleatoria proviene de una distribución normal. En este caso se cuenta con las siguientes hipótesis:

$H_0$: Los datos provienen de una distribución normal

$H_1$: Los datos no provienen de una distribución normal

El test provee como resultados un estadístico y un p-valor. Cuando el p-valor obtenido es superior a cierto umbral --usualmente 0.05--, decimos que no podemos rechazar la hipótesis de que los datos provienen de una distribución normal.

La implementación utilizada (scipy) indica que para muestras de más de 5000 números, se pierde precisión en los valores del test. Por lo que se realizará el test para muestras de tamaños N = 100, 1000 y 5000. (Ver: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html?highlight=shapiro)

In [31]:
m1_sw, m2_sw, m3_sw = [], [], []     #muestras
for i in range(100):
  z1 = generar_normal_0_1()
  z2 = transformar_gaussiana(z1, 10, 2)
  m1_sw.append(z2)

for i in range(1000):
  z1 = generar_normal_0_1()
  z2 = transformar_gaussiana(z1, 10, 2)
  m2_sw.append(z2)

for i in range(5000):
  z1 = generar_normal_0_1()
  z2 = transformar_gaussiana(z1, 10, 2)
  m3_sw.append(z2)
In [32]:
figsw, (ax1sw, ax2sw, ax3sw) = plt.subplots(1, 3, figsize=(20, 5))
figsw.tight_layout(pad=3)
ax1sw.hist(m1_sw, 10)
ax1sw.set_xlabel('Valores generados', fontsize=15)
ax1sw.set_title("Histograma (n = 100)", fontsize=20)
ax2sw.hist(m2_sw, 25)
ax2sw.set_xlabel('Valores generados', fontsize=15)
ax2sw.set_title("Histograma (n = 1000)", fontsize=20)
ax3sw.hist(m3_sw, 50)
ax3sw.set_xlabel('Valores generados', fontsize=15)
ax3sw.set_title("Histograma (n = 5000)", fontsize=20)

plt.show();
In [33]:
from scipy import stats
shapiro1 = stats.shapiro(m1_sw)
shapiro2 = stats.shapiro(m2_sw)
shapiro3 = stats.shapiro(m3_sw)

display ("Valor del estadístico del test Shapiro-Wilk con n = 100: " + str(shapiro1.statistic))
display ("P-valor del test Shapiro-Wilk con n = 100: " + str(shapiro1.pvalue))
display ("Valor del estadístico del test Shapiro-Wilk con n = 1000: " + str(shapiro2.statistic))
display ("P-valor del test Shapiro-Wilk con n = 1000: " + str(shapiro2.pvalue))
display ("Valor del estadístico del test Shapiro-Wilk con n = 5000: " + str(shapiro3.statistic))
display ("P-valor del test Shapiro-Wilk con n = 5000: " + str(shapiro3.pvalue))
'Valor del estadístico del test Shapiro-Wilk con n = 100: 0.9904636740684509'
'P-valor del test Shapiro-Wilk con n = 100: 0.7020654678344727'
'Valor del estadístico del test Shapiro-Wilk con n = 1000: 0.9989786744117737'
'P-valor del test Shapiro-Wilk con n = 1000: 0.8619867563247681'
'Valor del estadístico del test Shapiro-Wilk con n = 5000: 0.9995430111885071'
'P-valor del test Shapiro-Wilk con n = 5000: 0.28750666975975037'

Dado que el p-valor es superior a 0.05, no rechazamos la hipótesis de que los datos provienen de una distribución normal.

Proceso de Poisson¶

In [34]:
import os

wd = os.getcwd()
if 'TP1' not in wd: wd += '/TP1'

arribos_real = []
with open(f"{wd}/tiempos_entre_arribos.txt") as f:
    for l in f.readlines(): arribos_real.append(float(l.strip()))

tasa = len(arribos_real)/sum(arribos_real)
print(f"Tengo {len(arribos_real)} arribos en {sum(arribos_real)} horas")
print(f"Mi tasa de arribos empírica es de {tasa} arribos por hora")
Tengo 10000 arribos en 1011.1994386169154 horas
Mi tasa de arribos empírica es de 9.889245996494681 arribos por hora
In [35]:
# Para generar eventos de Poisson:
#  1. Generar muestras de una distribución exponencial Z, 
#     con parametro Lambda: z1, z2, z3...
#  2. Definir tiempos de eventos conteo de Poisson: 
#     t0 = 0, t1 = t0 + z1, t2 = t1 + z2....

def iter_exponencial(lmbda, generador, n):
    """Genera una secuencia que siga una distribución exponencial, 
       dado un generador que sigue una distribución uniforme"""
    u = np.array([next(generador) for i in range(int(n))])
    x = -np.log(1-u)/lmbda
    return iter(x)

def generar_poisson(limite, lmbda, generador):
    expo = iter_exponencial(lmbda, generador, lmbda*limite*1.5)
    arribos = []
    t = 0
    while t <= limite:
        arribo = next(expo)
        arribos.append(arribo)
        t += arribo
    return arribos

limite = 24 * 31 # Un mes
arribos = generar_poisson(limite, tasa, LXM(42))

print(f"Simulé {len(arribos)} arribos en {sum(arribos)} horas")
print(f"Mi tasa de arribos simulada es de {len(arribos)/sum(arribos)} arribos por hora")
Simulé 7300 arribos en 744.064522714941 horas
Mi tasa de arribos simulada es de 9.810977109033201 arribos por hora
In [36]:
def comparar_poissones(poissones_con_label, n = 15):
    for arribos, label in poissones_con_label:
        plt.step(np.cumsum(arribos[:n]),range(n),where= 'post',label=label)
    plt.legend()
    plt.title(f"Comparación de los primeros {n} arribos")
    plt.show()
    
comparar_poissones([(arribos_real, "Poisson Empirica"), (arribos, "Poisson Simulada")])
In [38]:
limite = 96
poissones = []
for i in range(1000):
    poissones.append(generar_poisson(limite, tasa, LXM(i)))
    
comparar_poissones([(p, f"Poisson Generada #{i}") for i, p in enumerate(poissones[:5])])
In [39]:
import math

# PMF de poisson, a mano, para comparar con los valores teóricos
pmf = lambda x, lmbda: ((lmbda ** x) * (math.e ** -lmbda)) / (math.factorial(x))

print("Probabilidad que el primer vehículo arribe antes de los 10 minutos.")
diez_min = 1/6


# La probabilidad de que haya al menos un arribo en 10 minutos equivale a
#   1 - la probabilidad de que no haya ningun arribo en 10 minutos
print(f"* Probabilidad Teórica => {1 - pmf(0, tasa*diez_min)}")

favorable = []
for p in poissones:
    favorable.append(p[0] <= diez_min)

print(f"* Probabilidad Simulada => {sum(favorable)/len(favorable)}")
Probabilidad que el primer vehículo arribe antes de los 10 minutos.
* Probabilidad Teórica => 0.8076055651533732
* Probabilidad Simulada => 0.818
In [40]:
print("Probabilidad que el undécimo vehículo arribe después de los 60 minutos.")
sesenta_min = 1

# La probabilidad de que el vehiculo 11 llegué despues de 60 minutos equivale a
#   la sumatoria de que hayan llegado 0,1,2,...,10 vehiculos en 60 minutos
# Ojo con los off-by-one! `for i in range(11)` === de 0 a 10
print(f"* Probabilidad Teórica => {sum([pmf(i, tasa*sesenta_min) for i in range(11)])}")

favorable = []
for p in poissones:
    arribos_acumulados = np.cumsum(p)
    # El arribo en el indice [10] === El undecimo arribo
    favorable.append(arribos_acumulados[10] >= 1)

print(f"* Probabilidad Simulada => {sum(favorable)/len(favorable)}")
Probabilidad que el undécimo vehículo arribe después de los 60 minutos.
* Probabilidad Teórica => 0.596893339441658
* Probabilidad Simulada => 0.579
In [41]:
# Como nuestra pmf a mano da overflow al calcular un número grande como lambda a la 750,
#   importamos directamente el modulo de scipy
from scipy.stats import poisson

print("Probabilidad que arriben al menos 750 vehículos antes de las 72 horas.")
setentidos_horas = 72

# La probabilidad de que arriben al menos 750 vehiculos en 72 horas equivale a
#   1 - la sumatoria de que lleguen 0,1,...,749 en 72 horas
print(f"* Probabilidad Teórica => {1 - sum([poisson.pmf(i, tasa*setentidos_horas) for i in range(750)])}")

favorable = []
for p in poissones:
    arribos_acumulados = np.cumsum(p)
    favorable.append(arribos_acumulados[749] <= 72)

print(f"* Probabilidad Simulada => {sum(favorable)/len(favorable)}")
Probabilidad que arriben al menos 750 vehículos antes de las 72 horas.
* Probabilidad Teórica => 0.08097848535892749
* Probabilidad Simulada => 0.084

Ejercicio 5¶

A partir del generador de número al azar implementado en el ejercicio 1, y del dataset provisto, obtenido del sitio de datos abiertos del Gobierno de la Ciudad de Buenos Aires (data.buenosaires.gob.ar), el cual contiene información geográfica de barrios de la Ciudad de Buenos Aires, se pide:

  • Proponer e implementar un método que permita generar coordenadas (latitud y longitud) distribuidas de forma uniforme en cada uno de los barrios.
  • Graficar los puntos generados en el mapa
In [44]:
class EJ5:
    
    def __init__(self, gdf, lxm):
        self.gdf = gdf.copy()
        self.lxm = lxm
        self.puntos_en_barrios = {}
        for barrio in self.gdf['BARRIO'].tolist():
            self.puntos_en_barrios[barrio] = [[],[]]

    def emitir_barrios(self):
        display(self.gdf['BARRIO'])

    #metodo auxiliar para trasladar un barrio al primer cuadrante
    def trasladar_barrio(self, datos_barrio):
        extremos_barrio = datos_barrio.geometry.bounds

        signo1 = 1
        signo2 = 1

        if(extremos_barrio.iloc[0]['minx'] < 0):
            signo1 = -1
        if(extremos_barrio.iloc[0]['miny'] < 0):
            signo2 = -1

        return datos_barrio.translate(extremos_barrio.iloc[0]['minx']*signo1, extremos_barrio.iloc[0]['miny']*signo2)

    #Generar n puntos en un barrio
    def generar_en_barrio(self, barrio, cantidad):
        x, y = [], []

        #obtengo la geometria del barrio, la traslado al primer cuadrante
        #y busco el maximo entre sus extremos para formar un cuadrado que los contenga
        datos_barrio = self.gdf[self.gdf['BARRIO']==barrio]
        g_barrio = datos_barrio['geometry']
        g_barrio_t = self.trasladar_barrio(datos_barrio)
        extremos_b = datos_barrio.bounds.to_numpy()
        extremos_t = g_barrio_t.bounds.to_numpy()
        maximo_t = np.amax(extremos_t)

        #genero puntos dentro del cuadrado que obtuve, y los acepto si están dentro del barrio
        x = []
        y = []
        while len(x) < cantidad:
            z1 = next(self.lxm) * maximo_t
            z2 = next(self.lxm) * maximo_t
            ps = shapely.geometry.Point(z1, z2)
 
            if g_barrio_t.contains(ps).any():
                x.append(z1 + datos_barrio.geometry.bounds.minx)
                y.append(z2 + datos_barrio.geometry.bounds.miny)

        self.puntos_en_barrios[barrio] = [x, y]


    #Generar n puntos en todos los barrios
    def generar_en_todos(self, n):
        for barrio in self.gdf['BARRIO'].tolist():
            self.generar_en_barrio(barrio, n)

    #Graficar un barrio y los puntos generados en el mismo
    def graficar_en_barrio(self, barrio):
        x = self.puntos_en_barrios[barrio][0]
        y = self.puntos_en_barrios[barrio][1]
        title_string = "Puntos generados en el barrio de {} (n = {})".format(barrio, len(x))
        
        fig, ax = plt.subplots(figsize=(10, 10))
        fig.tight_layout(pad=3)
        ax.set_title(title_string, fontsize=20)

        datos_barrio = self.gdf.loc[self.gdf['BARRIO']==barrio]
        datos_barrio.plot(ax=ax, color='white', edgecolor='black')
        ax.scatter(x, y, s=3, color="red")
        ax.grid()

        return fig, ax

    #Graficar el mapa completo y los puntos generados en un barrio
    def graficar_mapa_completo_barrio(self, nombre_barrio):
        x = self.puntos_en_barrios[nombre_barrio][0]
        y = self.puntos_en_barrios[nombre_barrio][1]
        title_string = "Puntos generados en el barrio de {} (n = {})".format(nombre_barrio, len(x))
        
        fig, ax = plt.subplots(figsize=(10, 10))
        fig.tight_layout(pad=3)
        ax.set_title(title_string, fontsize=20)
        self.gdf.plot(ax=ax, color='white', edgecolor='black')
        ax.scatter(x, y, s=3, color="red")
        ax.grid()
        
        return fig, ax

    #Graficar el mapa completo y los puntos generados en todos los barrios (si hay)
    def graficar_mapa_completo_todos(self):
        title_string = "Puntos generados en todos los barrios."
        fig, ax = plt.subplots(figsize=(10, 10))
        fig.tight_layout(pad=3)
        ax.set_title(title_string, fontsize=20)

        self.gdf.plot(ax=ax, color='white', edgecolor='black')

        for b in self.gdf['BARRIO'].tolist():
            x = self.puntos_en_barrios[b][0]
            y = self.puntos_en_barrios[b][1]
            ax.scatter(x, y, s=3, color="red")
        ax.grid()

        return fig, ax

    def help(self):
        display("generar_en_barrio('barrio', cantidad)")
        display("graficar_mapa_completo_barrio('barrio')")
        display("graficar_mapa_completo_todos()")
        display("graficar_en_barrio('barrio')")
        display("generar_en_todos(cantidad)")
  
In [45]:
ej5test = EJ5(geodataframe, LXM(42))
ej5test.generar_en_barrio('PUERTO MADERO', 100)
ej5test.generar_en_todos(100)
In [46]:
ej5test.graficar_en_barrio('PUERTO MADERO')
Out[46]:
(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)
In [47]:
ej5test.graficar_mapa_completo_barrio('PUERTO MADERO')
Out[47]:
(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)
In [48]:
ej5test.graficar_mapa_completo_todos()
Out[48]:
(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en todos los barrios.'}>)

Comparaciones con un GCL¶

In [49]:
# Ahora, rehagamos partes del trabajo, pero basados en un GCL en vez de en LXM
# Son los resultados distintos?
padrones = [97568, 100029, 100033, 104030]
semilla = int(np.mean(padrones))

def GCL(x, a = 1013904223, c = 1664525, m = 232):
    while True:
        x = ((a * x + c) % m) / m
        yield x
In [50]:
n = 5
g = GCL(semilla) 
print(f"{n} números generados al azar, en el intervalo 0-1:")
for r, i in zip(g,range(n)):
    print(f"* {r}")
5 números generados al azar, en el intervalo 0-1:
* 0.9224137931034483
* 0.29819411399035617
* 0.5019334719098848
* 0.895838293535956
* 0.008367325211393422
In [52]:
scatterplot.show()
In [53]:
Chi2(numero_de_muestras=10000, g=GCL(semilla)).realizar_prueba()
Frecuencias: [1005, 996, 1013, 1018, 1002, 992, 987, 1020, 994, 973]
Estadístico: 1.94
Límite superior: 11.07
El test acepta la hipotesis nula
In [54]:
limite = 24 * 31 # Un mes
arribos_gcl = generar_poisson(limite, tasa, GCL(semilla))

print(f"Simulé {len(arribos_gcl)} arribos en {sum(arribos_gcl)} horas")
print(f"Mi tasa de arribos simulada es de {len(arribos_gcl)/sum(arribos_gcl)} arribos por hora")
comparar_poissones([(arribos_real, "Poisson Empirica"), (arribos, "Poisson Simulada con LXM"), (arribos_gcl, "Poisson Simulada con GCL")])
Simulé 7394 arribos en 744.0969199691178 horas
Mi tasa de arribos simulada es de 9.936877578134409 arribos por hora
In [55]:
ej5test = EJ5(geodataframe, GCL(semilla))
ej5test.generar_en_barrio('PUERTO MADERO', 100)
ej5test.generar_en_todos(100)
In [56]:
ej5test.graficar_en_barrio('PUERTO MADERO')
Out[56]:
(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)
In [57]:
ej5test.graficar_mapa_completo_barrio('PUERTO MADERO')
Out[57]:
(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)
In [58]:
ej5test.graficar_mapa_completo_todos()
Out[58]:
(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en todos los barrios.'}>)

Conclusión: Los resultados con LXM y con el GCL a lo largo del TP nos parecen bastante similares, lo cual nos parece que tiene muchísimo sentido! LXM se basa fuertemente en un GCL.

Trabajo Práctico 2¶

Page Rank de Actores y Actrices¶

In [59]:
import os
from bs4 import BeautifulSoup

wd = os.getcwd()
if 'TP2' not in wd: wd += '/TP2'

DIR = f'{wd}/paginas'
PAGINAS = os.listdir(DIR)

links = [] # Lista de tuplas de src->dst
palabras = {} # Dict de pagina: palabras

for f in PAGINAS:
    with open(os.path.join(DIR, f)) as file:
        soup = BeautifulSoup(file.read())
        palabras[f] = soup.get_text().split()

        for link in soup.find_all('a'):
            dst = link.get('href').replace('http://','')
            links.append((f,dst))

links
Out[59]:
[('jonvoight.html', 'angelinajolie.html'),
 ('jonvoight.html', 'angelinajolie.html'),
 ('jonvoight.html', 'bradpitt.html'),
 ('bradpitt.html', 'jenniferaniston.html'),
 ('bradpitt.html', 'angelinajolie.html'),
 ('bradpitt.html', 'martinscorcese.html'),
 ('bradpitt.html', 'angelinajolie.html'),
 ('robertdeniro.html', 'martinscorcese.html'),
 ('angelinajolie.html', 'jonvoight.html'),
 ('angelinajolie.html', 'bradpitt.html')]
In [ ]:
 
In [ ]:
 
In [61]:
import networkx as nx

G = nx.DiGraph()
G.add_edges_from(links)
nx.draw_shell(G, with_labels=True, node_color='white', edge_color='grey')
In [62]:
matriz_ady = nx.to_pandas_adjacency(G)
matriz_ady
Out[62]:
jonvoight.html angelinajolie.html bradpitt.html jenniferaniston.html martinscorcese.html robertdeniro.html
jonvoight.html 0.0 1.0 1.0 0.0 0.0 0.0
angelinajolie.html 1.0 0.0 1.0 0.0 0.0 0.0
bradpitt.html 0.0 1.0 0.0 1.0 1.0 0.0
jenniferaniston.html 0.0 0.0 0.0 0.0 0.0 0.0
martinscorcese.html 0.0 0.0 0.0 0.0 0.0 0.0
robertdeniro.html 0.0 0.0 0.0 0.0 1.0 0.0
In [63]:
suma_filas = matriz_ady.sum(axis=1)
matriz_page_rank = matriz_ady.div(suma_filas, axis=0)
matriz_page_rank.fillna(1/len(PAGINAS), inplace=True)
matriz_page_rank
Out[63]:
jonvoight.html angelinajolie.html bradpitt.html jenniferaniston.html martinscorcese.html robertdeniro.html
jonvoight.html 0.000000 0.500000 0.500000 0.000000 0.000000 0.000000
angelinajolie.html 0.500000 0.000000 0.500000 0.000000 0.000000 0.000000
bradpitt.html 0.000000 0.333333 0.000000 0.333333 0.333333 0.000000
jenniferaniston.html 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667
martinscorcese.html 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667
robertdeniro.html 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000
In [64]:
import numpy as np

# Hacemos N iteraciones de la matriz
iteraciones = np.linalg.matrix_power(matriz_page_rank, 50)

page_rank = {p: iteraciones[0][i] for i, p in enumerate(matriz_page_rank.index)}
page_rank
Out[64]:
{'jonvoight.html': 0.16216216216216203,
 'angelinajolie.html': 0.21621621621621548,
 'bradpitt.html': 0.2432432432432427,
 'jenniferaniston.html': 0.13513513513513484,
 'martinscorcese.html': 0.18918918918918876,
 'robertdeniro.html': 0.05405405405405395}
In [65]:
import random

BUSQUEDAS = 20

todas_palabras = set()
for p in palabras.values():
    todas_palabras.update(p)

for p in random.choices(list(todas_palabras), k=BUSQUEDAS):
    matches = [f for f in PAGINAS if p in palabras[f]]
    print(f"'{p}' aparece en: {sorted(matches, key=lambda x: -page_rank[x])}")
'Rumor' aparece en: ['jenniferaniston.html']
'Interrupted' aparece en: ['angelinajolie.html']
'recent' aparece en: ['jonvoight.html']
'man' aparece en: ['jonvoight.html']
'Four' aparece en: ['bradpitt.html']
'1994' aparece en: ['bradpitt.html']
'Twelve' aparece en: ['bradpitt.html']
'has' aparece en: ['bradpitt.html', 'angelinajolie.html', 'martinscorcese.html', 'jonvoight.html', 'jenniferaniston.html']
'Zahara,' aparece en: ['bradpitt.html', 'angelinajolie.html']
'end' aparece en: ['jonvoight.html']
'Scorsese' aparece en: ['martinscorcese.html']
'Midnight' aparece en: ['jonvoight.html']
'portrayal' aparece en: ['angelinajolie.html']
'United' aparece en: ['bradpitt.html']
'World' aparece en: ['martinscorcese.html']
'She' aparece en: ['angelinajolie.html', 'jenniferaniston.html']
'Interview' aparece en: ['bradpitt.html']
'Maddox,' aparece en: ['bradpitt.html', 'angelinajolie.html']
'an' aparece en: ['bradpitt.html', 'angelinajolie.html', 'martinscorcese.html', 'jonvoight.html', 'jenniferaniston.html', 'robertdeniro.html']
'Smith' aparece en: ['bradpitt.html', 'angelinajolie.html']

Ejercicio 2¶

En base a lo presentado en el trabajo “Queuing theory application in imaging service analysis for small Earth observation satellites”, simular los resultados obtenidos sobre la longitud de imágenes en cola esperando ser procesadas en la sección 3.1. Pure image capture service system El ejercicio se puede resolver utilizando simpy o programación tradicional (a elección del grupo)

Notas sobre el paper¶

El paper trata de la aplicación de la teoría de colas al procesamiento de solicitudes de captura de imágenes satelitales. En particular el modelo descrito consta de dos etapas: la captura de imágenes y la descarga de imágenes. A continuación intentaremos replicar los resultados obtenidos al modelar el servicio que procesa solicitudes de captura de imágenes, descrito en la sección 3.1 del paper.

El servicio de captura de imágenes inicia con al llegada de una solicitud, y termina con la captura de la imágen solicitada.

Las solicitudes de captura de imágenes ingresan a la cola de acuerdo con un proceso de poisson de tasa $\lambda$.

Adicionalmente el modelo supone lo siguiente:

  • No hay límite a la longitud de la cola de solicitudes.
  • La politica de atención de solicitudes es FOFS (se priorizan las imagenes de las ubicaciones más cercanas en la ruta del satélite).
  • Las ubicaciones de las solicitudes están distribuidas de manera uniforme sobre el planeta.
In [67]:
class Metricas():
    def __init__(self, tasaArribos, tasaServicio):
        self.tasaArribos = tasaArribos
        self.tasaServicio = tasaServicio
        self.arribos = []
        self.salidas = []
        self.tiempoEnCola = []
        self.tiemposTotales = []
        self.cantidadEnCola = []
        self.tiemposDeServicio = []
        #self.tiemposDeServicioCola = []
        self.tasasServicio = []
        self.longitudMediaTeorica = 0
        self.longitudMediaEmpirica = 0
        self.tiempoMedioTeorico = 0
        self.tiempoMedioEmpirico = 0
        self.tasaServicioMedia = 0
        self.probabilidades = []

    def calcular_probabilidades(self):
        for i in range(0, np.amax(self.cantidadEnCola)):
            p = (self.tasaArribos/self.tasaServicio)**i
            p = p*np.e**-(self.tasaArribos/self.tasaServicio)
            p = p/np.math.factorial(i)
            self.probabilidades.append(p)

    def calcular_metricas_paper(self):
        self.calcular_probabilidades()

        #longitud media teorica
        for i in range(len(self.probabilidades)):
            self.longitudMediaTeorica = self.longitudMediaTeorica + i*self.probabilidades[i]
        
        #tiempo de espera medio teorico
        self.tiempoMedioTeorico = self.longitudMediaTeorica / self.tasaArribos

        #valores calculados a partir de la simulacion
        self.longitudMediaEmpirica = np.mean(self.cantidadEnCola)
        self.tiempoMedioEmpirico = np.mean(self.tiemposDeServicio)
        self.tasaServicioMedia = np.mean(self.tasasServicio)
In [68]:
#generador para modelar la llegada de solicitudes
def solicitud(env, servicio, cant_arribos, tasaArribos, tasaServicio, metricas):
    solID = 0
    while solID < cant_arribos:
        tiempoSolicitud = generar_exponencial(tasaArribos)
        yield env.timeout(tiempoSolicitud)
        
        solID += 1
        metricas.arribos.append(env.now)
        env.process(captura_imagen(env, servicio, solID, tasaServicio, metricas))
In [69]:
#generador para modelar el servicio de captura de imágenes
def captura_imagen(env, servicio, solID, tasaServicio, metricas):
    with servicio.request() as req:
        inicio_espera = env.now #esto estaba antes de yield req
        yield req
        fin_espera = env.now

        n = max(1,len(servicio.queue)) #esto estaba despues de yield req
        tiempoDeServicio = generar_exponencial(tasaServicio*n)
        metricas.tiemposDeServicio.append(tiempoDeServicio)
        metricas.cantidadEnCola.append(len(servicio.queue))
        metricas.tasasServicio.append(tasaServicio*n)

        inicio_servicio = env.now
        yield env.timeout(tiempoDeServicio)
        fin_servicio = env.now

        metricas.tiempoEnCola.append(fin_espera - inicio_espera)
        metricas.tiemposTotales.append(fin_servicio - inicio_espera)
        metricas.salidas.append(env.now)
In [70]:
#inicio la simulación
def simular(cant_satelites, cant_arribos, tasaArribos, tasaServicio, semilla):
    print(f'Simulando {cant_arribos} arribos con lambdaArribos = {tasaArribos} y lambdaServicio = {tasaServicio}')
  
    np.random.seed(seed=semilla)
    metricas = Metricas(tasaArribos, tasaServicio)

    env = simpy.Environment()
    servicio = simpy.Resource(env, capacity = cant_satelites)
    env.process(solicitud(env, servicio, cant_arribos, tasaArribos, tasaServicio, metricas))
    env.run()
    return metricas
In [71]:
metricasHist = []
a_arribos = np.linspace(1, 15, 15)
a_llegadas = np.array([1])
longitudes = []
tiempos = []
longitudMediaTasaMedia = []
for i in a_arribos:
    for j in a_llegadas:
        metricas = simular(cant_satelites=1, cant_arribos=200, tasaArribos=i, tasaServicio=j, semilla=2951123)
        metricas.calcular_metricas_paper()
        longitudes.append((metricas.longitudMediaTeorica, metricas.longitudMediaEmpirica))
        tiempos.append((metricas.tiempoMedioTeorico, metricas.tiempoMedioEmpirico))
        longitudMediaTasaMedia.append((metricas.longitudMediaEmpirica, metricas.tasaServicioMedia))
        metricasHist.append(metricas)
Simulando 200 arribos con lambdaArribos = 1.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 2.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 3.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 4.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 5.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 6.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 7.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 8.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 9.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 10.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 11.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 12.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 13.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 14.0 y lambdaServicio = 1
Simulando 200 arribos con lambdaArribos = 15.0 y lambdaServicio = 1

1.) Verificaciones sobre los procesos de Poisson

1.1.) Evolución de arribos y salidas.

Tomando tres simulaciones de ejemplo, con distintos valores de tasas de arribo de solicitudes (λ=1.0, λ=5.0, λ=9.0), graficamos la evolución de los arribos y salidas de solicitudes.

En particular notamos como la tasa de atención del servicio de captura de imágenes se ajusta en función de la cantidad de solicitudes en cola, de forma que, a pesar de que la tasa de arribos es superior a la tasa de atención, los arribos y salidas evolucionan a la par. En condiciones normales esperaríamos que el gráfico de arribos tenga una pendiente mucho mayor a la del gráfico de salidas.

In [72]:
def graficar_arribos_salidas(m1, m2, m3):
    fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))
    ax[0].step(x=m1.arribos, y=range(len(m1.arribos)))
    ax[0].step(x=m1.salidas, y=range(len(m1.salidas)))
    ax[1].step(x=m2.arribos, y=range(len(m2.arribos)))
    ax[1].step(x=m2.salidas, y=range(len(m2.salidas)))
    ax[2].step(x=m3.arribos, y=range(len(m3.arribos)))
    ax[2].step(x=m3.salidas, y=range(len(m3.salidas)))
    fig.suptitle("Evolución de arribos y salidas", fontsize=20)
    ax[0].set_xlabel(f'Arribos y salidas \n (λ = {m1.tasaArribos}, μ = {m1.tasaServicio})', fontsize=15)
    ax[1].set_xlabel(f'Arribos y salidas \n (λ = {m2.tasaArribos}, μ = {m2.tasaServicio})', fontsize=15)
    ax[2].set_xlabel(f'Arribos y salidas \n (λ = {m3.tasaArribos}, μ = {m3.tasaServicio})', fontsize=15)
    ax[0].legend([f'λ = {m1.tasaArribos}, μ = {m1.tasaServicio}'])
    ax[1].legend([f'λ = {m2.tasaArribos}, μ = {m2.tasaServicio}'])
    ax[2].legend([f'λ = {m3.tasaArribos}, μ = {m3.tasaServicio}'])
    ax[2].grid()
    ax[1].grid()
    ax[0].grid()
    plt.show()
    
graficar_arribos_salidas(metricasHist[0], metricasHist[4], metricasHist[8])

1.2. Densidad de los tiempos de servicio

Para las simulaciones seleccionadas en el punto anterior, graficamos la densidad de los tiempos de servicio. Verificamos que, para valores de λ muy superiores a la tasa de atención, los cuales se corresponden con una mayor cantidad de solicitudes en cola, los tiempos de servicio se acumulan en torno a 0.0, es decir, son menores que para valores de λ cercanos a la tasa de atención.

In [73]:
def graficar_densidad_tiempos_servicio(m1, m2, m3):
    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(15,5))
    p = sns.kdeplot(m1.tiemposDeServicio,ax=ax)
    p = sns.kdeplot(m2.tiemposDeServicio,ax=ax)
    p = sns.kdeplot(m3.tiemposDeServicio,ax=ax)
    ax.legend([f'λ = {m1.tasaArribos}, μ = {m1.tasaServicio}', f'λ = {m2.tasaArribos}, μ = {m2.tasaServicio}', f'λ = {m3.tasaArribos}, μ = {m3.tasaServicio}'])
    ax.set_xticks(np.arange(0, 7, 0.2))
    ax.set_yticks(np.arange(0, 4.6, 0.2))
    ax.grid()
    #fig.tight_layout()
    fig.suptitle("Densidad de tiempos de servicio", fontsize=20)

graficar_densidad_tiempos_servicio(metricasHist[0], metricasHist[4], metricasHist[8])

2.) Resultados del modelo

2.1.) Relación entre la tasa de atención y la longitud de la cola.

Debido a la politica de priorización de solicitudes, el modelo predice que habrá una relación lineal entre la tasa de atención de solicitudes y el largo de la cola. A continuación graficamos la correspondencia entre ambas cantidades, para todas las simulaciones realizadas:

In [74]:
fig, ax5 = plt.subplots(nrows=1, ncols=1, figsize=(10, 10))
ax5.scatter(*zip(*longitudMediaTasaMedia))
ax5.set_xlabel('Longitud media de la cola de solicitudes', fontsize=15)
ax5.set_ylabel('Tasa media del servicio de captura de imágenes  (sol./tiempo)', fontsize=15)
ax5.set_title("Longitudes de cola medias", fontsize=20)
ax5.set_xticks(np.arange(0,30, 1))
ax5.set_yticks(np.arange(0,30, 1))
ax5.grid()

2.2) Probabilidades de las posibles longitudes de la cola

El modelo da una expresión para las probabilidades estacionarias de las longitudes de cola.

$$p_n = \frac{(λ/μ_0)e^{-(λ/μ_0)}}{n!} (n≥0)$$

Analizamos para los casos vistos anteriormente (λ=1.0, λ=5.0, λ=9.0) si la densidad de cantidad de solicitudes en cola obtenida de forma empírica se asemeja a los valores dados por el modelo.

In [75]:
def grafico_probabilidades(m1, m2, m3):
    met = [m1, m2, m3]
    fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(25, 5))
    for i in range(0, 3):
        p = sns.kdeplot(met[i].cantidadEnCola,ax=ax[i], c="red")
        ax[i].scatter(y=met[i].probabilidades, x=range(len(met[i].probabilidades)))
        ax[i].legend([f'λ = {met[i].tasaArribos}, μ = {met[i].tasaServicio}', 'Probabilidades Teóricas'])
        ax[i].set_xlabel('Solicitudes en cola', fontsize=15)
        ax[i].set_ylabel('Densidad', fontsize=15)
        ax[i].set_xticks(np.arange(0,25, 1))
        ax[i].grid()
    fig.suptitle("Densidad de cantidad de solicitudes en cola", fontsize=20)

grafico_probabilidades(metricasHist[0], metricasHist[4], metricasHist[8])

2.3) Longitud media de la cola y tiempo medio de espera

El tercer resultado es que el modelo da expresiones para la longitud media de la cola ($L_i$), y para el tiempo medio que una solicitud espera para ser atendida ($W$).

$$L_i = ∑{n p_n}$$$$$$$$W = \frac{L_i}{λ} $$

Graficamos los valores teóricos y los obtenidos a partir de las simulaciones, esperando obtener una relacion lineal entre ellos:

In [76]:
fig, ax6 = plt.subplots(nrows=1, ncols=2, figsize=(20, 10))
ax6[0].scatter(*zip(*longitudes))
ax6[0].set_xlabel('Longitud Media Teórica', fontsize=15)
ax6[0].set_ylabel('Longitud Media Empirica', fontsize=15)
ax6[0].set_title("Longitudes de cola medias", fontsize=20)
ax6[0].grid()
ax6[0].set_xticks(np.arange(0,20, 1))
ax6[0].set_yticks(np.arange(0,20, 1))
ax6[1].scatter(*zip(*tiempos))
ax6[1].set_xlabel('Tiempo medio de servicio Teórico', fontsize=15)
ax6[1].set_ylabel('Tiempo medio de servicio Empírico', fontsize=15)
ax6[1].set_title("Tiempos medios de servicio", fontsize=20)
ax6[1].grid()
ax6[1].set_xticks(np.arange(0.9,1.1, 0.01))
ax6[1].set_yticks(np.arange(0,0.3, 0.01))
plt.show()

Ejercicio 3: Diseño del Web Service¶

In [78]:
class Consulta:
    def __init__(self, tiempo_arribo, duracion, tiempo_finalizacion_anterior):
        if tiempo_arribo < tiempo_finalizacion_anterior:
            self.duracion_espera = tiempo_finalizacion_anterior - tiempo_arribo
            self.tuvo_que_esperar = True
        else:
            self.duracion_espera = 0
            self.tuvo_que_esperar = False
                            
        self.tiempo_finalizacion = tiempo_arribo + duracion + self.duracion_espera
In [79]:
def simular_consultas(tiempos_de_arribo, duraciones):
    consultas = []

    # Uso -1 como tiempo de finalizacion del anterior para que empieze ni bien llega
    primera_consulta = Consulta(tiempos_de_arribo[0], duraciones[0], tiempo_finalizacion_anterior = -1)
    consultas.append(primera_consulta)

    for tiempo_arribo, duracion in zip(tiempos_de_arribo[1:], duraciones[1:]):
        consulta = Consulta(tiempo_arribo,
                        duracion, 
                        tiempo_finalizacion_anterior = consultas[-1].tiempo_finalizacion)
        consultas.append(consulta)
        
    tiempo_espera_total = 0
    cant_consultas_que_no_esperaron = 0
    for consulta in consultas:
        tiempo_espera_total += consulta.duracion_espera
        if not consulta.tuvo_que_esperar:
            cant_consultas_que_no_esperaron += 1

    return tiempo_espera_total, cant_consultas_que_no_esperaron

Centralizado¶

In [82]:
t, f = simular_consultas_centralizado(100000, 4, 0.8)
mostrar_resultados("Base de datos central", t,f, 100000)
--------------Base de datos central--------------
	El tiempo de espera promedio fue 11 Horas 3 Minutos 16.0 Segundos
	El 2e-05% de las consultas no tuvo que esperar

Distribuido¶

In [84]:
t1, f1, l1, t2, f2, l2 = simular_consultas_distribuido(100000, 4, 0.7, 1, 0.7)
mostrar_resultados("Ambas bases distribuidas", t1 + t2, f2 + f2, l1 + l2)
mostrar_resultados("Base de datos distribuida 1", t1, f1, l1)
mostrar_resultados("Base de datos distribuida 2", t2, f2, l2)
--------------Ambas bases distribuidas--------------
	El tiempo de espera promedio fue 7 Horas 40 Minutos 21.0 Segundos
	El 4e-05% de las consultas no tuvo que esperar
--------------Base de datos distribuida 1--------------
	El tiempo de espera promedio fue 2 Horas 51 Minutos 34.0 Segundos
	El 7e-05% de las consultas no tuvo que esperar
--------------Base de datos distribuida 2--------------
	El tiempo de espera promedio fue 9 Horas 42 Minutos 30.0 Segundos
	El 3e-05% de las consultas no tuvo que esperar

Modificando los parametros¶

In [85]:
def simular_variaciones(cantidades_consultas=[100000], 
                        medias_tiempos_de_arribo=[4],
                        medias_duraciones_consultas = [0.8],  
                        medias_duraciones_consultas1 = [0.7],
                        medias_duraciones_consultas2 = [1],
                        lista_p = [0.7]):

    tiempos_centralizado = []
    frecuencias_centralizado = []

    tiempos_distribuido = []
    frecuencias_distribuido = []
                         
    for cant_consultas in cantidades_consultas:
        for media_tiempos_de_arribo in medias_tiempos_de_arribo:
            for media_duraciones_consultas in medias_duraciones_consultas:
                for media_duraciones_consultas1 in medias_duraciones_consultas1:
                     for media_duraciones_consultas2 in medias_duraciones_consultas2:
                            for p in lista_p:
                                t, f = simular_consultas_centralizado(cant_consultas, media_tiempos_de_arribo, 
                                                                      media_duraciones_consultas)
                                tiempos_centralizado.append(t / cant_consultas)
                                frecuencias_centralizado.append(f /cant_consultas)
    
                                t1, f1, l1, t2, f2, l2  = simular_consultas_distribuido(cant_consultas, 
                                                                                        media_tiempos_de_arribo,
                                                                                        media_duraciones_consultas1,
                                                                                        media_duraciones_consultas2,
                                                                                        p)
    
                                tiempos_distribuido.append((t1 + t2) / (l1 + l2))
                                frecuencias_distribuido.append((f1 + f2) / (l1 + l2))
            
            
    return tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido 

Media del tiempo de arribo¶

In [89]:
graficar_tiempo_espera(medias_tiempos_de_arribo, "Media de los tiempos de arribo de consultas",
                      tiempos_centralizado, tiempos_distribuido)
In [90]:
graficar_fraccion_que_no_espero(medias_tiempos_de_arribo, "Media de los tiempos de arribo de consultas",
                      frecuencias_centralizado, frecuencias_distribuido)

Lambda¶

In [92]:
graficar_tiempo_espera(lambda_l, "Lambda",
                      tiempos_centralizado, tiempos_distribuido, distribuido_es_referencia = True)

Lambda1¶

In [95]:
graficar_tiempo_espera(lambda_1, "Lambda1",
                      tiempos_centralizado, tiempos_distribuido, centralizado_es_referencia = True)

Lambda2¶

In [98]:
graficar_tiempo_espera(lambda_2, "Lambda2",
                      tiempos_centralizado, tiempos_distribuido, centralizado_es_referencia = True)

p¶

In [101]:
graficar_tiempo_espera(lista_p, "p",
                      tiempos_centralizado, tiempos_distribuido, centralizado_es_referencia = True)

Ejercicio 4¶

a. Simular 1000 días completos de 24 hrs.

b. Para un día en particular graficar la cantidad de billetes en el cajero luego de cada transacción.

c. Calcular el tiempo medio que los clientes demoraron en el sistema (espera + utilización del cajero)

d. ¿Recomienda a la entidad que implemente el cambio de cajero?

In [103]:
from scipy.stats import expon, uniform
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
from math import floor
from random import randint
import random
import matplotlib.pyplot as plt
from statistics import median


class Cajero:
    def __init__(self):
        self._billetes_en_cajero = 2000 # son todos de 100,  arranca el  día lleno
        self._capacidad_maxima = 2000  #billetes max q entran
        self._historial_dia_actual = []
        self._historial = []
        
    def nuevo_dia(self):
        self._billetes_en_cajero = 2000
        self._historial.append(self._historial_dia_actual)
        self._historial_dia_actual = []
    
    def retirar(self, cantidad_de_billetes):
        if self._billetes_en_cajero < cantidad_de_billetes:
            return False
        self._billetes_en_cajero -= cantidad_de_billetes
        self._historial_dia_actual.append(self._billetes_en_cajero)
        return True
        
    def depositar(self, cantidad_de_billetes):
        if self._billetes_en_cajero + cantidad_de_billetes > self._capacidad_maxima:
            return False
        self._billetes_en_cajero += cantidad_de_billetes
        self._historial_dia_actual.append(self._billetes_en_cajero)
        return True
        
    def historial(self):
        return self._historial
        
        
class Banco:
    def __init__(self):
        self._cajero = Cajero()
        self._minutos = 0 # al pasa 24 * 60 termina el día y restarteamos
        self._tiempos_de_espera = []
        self._clientes_perdidos = 0
        self._cantidad_de_clientes = 0
        self._cantidad_de_dias = 0
        # general
        self._minutos_en_hora = 24 *  60
        
    def _simular_dia(self):
        print(expon.rvs(size=3, random_state=10))
        # tarda en llegar n 
        
    def _grupo_cliente(self):
        return random.choices([1,2], weights=(75, 25))[0]
    
    def _resetear_dia(self):
        if self._minutos >= self._minutos_en_hora:
            self._minutos = self._minutos_en_hora - self._minutos
            self._cajero.nuevo_dia()
            self._cantidad_de_dias += 1
    
    def historial_de_cajero(self):
        return self._cajero.historial()
    
    def imprimir_estadisticas_de_simulacion(self):
        print(f'Pasaron {self._cantidad_de_dias} días')
        print(f'Cantidad de clientes: {self._cantidad_de_clientes}')
        print(f'Se perdieron {self._clientes_perdidos} clientes')
        porcentaje_clientes_perdidos = self._clientes_perdidos * 100 / self._cantidad_de_clientes
        print(f'Porcentaje de clientes perdidos: {floor(porcentaje_clientes_perdidos)}%')
        
    def imprimir_mediana(self):
        mediana = round(median(self._tiempos_de_espera), 2)
        print(f'Tiempo medio de espera: {mediana} minutos')
            
    def simular(self, dias, imprimir_estadisticas_clientes=False, imprimir_media_espera=False):
        tiempos_de_arribo = expon.rvs(size=144 * 60 * dias, scale=10, random_state=10)
        duraciones_1 = expon.rvs(size=len(tiempos_de_arribo), scale=1.5, random_state=10)
        duraciones_2 = expon.rvs(size=len(tiempos_de_arribo), scale=5, random_state=10)
        for i, tiempo in enumerate(tiempos_de_arribo):
            grupo_cliente =  self._grupo_cliente()
            self._minutos += tiempo
            duracion_de_persona = 0
            # Si es necesario se resetea el dia
            self._resetear_dia()
            if self._cantidad_de_dias == dias:
                break
            # Acciones de usuario
            if grupo_cliente == 1: # Retiran
                cant_billetes = randint(3, 50)
                fue_exitoso = self._cajero.retirar(cant_billetes)
                self._clientes_perdidos += 0 if fue_exitoso else 1
                duracion_de_persona = duraciones_1[i]
            else: # Depositan
                cant_billetes = randint(10, 110)
                fue_exitoso = self._cajero.depositar(cant_billetes)
                self._clientes_perdidos += 0 if fue_exitoso else 1
                duracion_de_persona = duraciones_2[i]
            demora = 0
            if i > 0 and tiempos_de_arribo[i-1] < duracion_de_persona:
                demora = duracion_de_persona - tiempos_de_arribo[i+1]
                self._minutos += demora
            self._tiempos_de_espera.append(duracion_de_persona + demora)
            self._cantidad_de_clientes += 1
        if imprimir_estadisticas_clientes:
            self.imprimir_estadisticas_de_simulacion()
        if imprimir_media_espera:
            self.imprimir_mediana()
    
In [104]:
# Simular 1000 días
banco = Banco()
banco.simular(dias=1000, imprimir_estadisticas_clientes=True)
    
Pasaron 1000 días
Cantidad de clientes: 159483
Se perdieron 2505 clientes
Porcentaje de clientes perdidos: 1%
In [105]:
# cantidad de billetes en el cajero dps de cada transacción
banco = Banco()
banco.simular(dias=1000)
monto_en_cajero = banco.historial_de_cajero()[1]

# make up some data
x = [i for i,monto in enumerate(monto_en_cajero)]

# plot
plt.scatter(x,monto_en_cajero, linewidths = 0.1)
# beautify the x-labels
plt.gcf().autofmt_xdate()
print(f"Cantidad de transacciones: {len(monto_en_cajero)}")
Cantidad de transacciones: 165
In [106]:
#  Calcular el tiempo medio que los clientes demoraron en el sistema (espera + utilización del cajero)
banco = Banco()
banco.simular(dias=1000, imprimir_media_espera=True)
Tiempo medio de espera: 1.05 minutos

¿Recomienda a la entidad que implemente el cambio de cajero?¶

Sí, ya que el promedio de clientes perdidos es entre 1% y 2%, mientras que nos indicaron antes que la perdida de clientes era de una 20%